{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TESTING = False # Use 1% of data for testing. IF false - full dataset takes more time\n",
    "MAX_FEATURE_SET = True\n",
    "#DIAGNOSIS = False\n",
    "IMPUTE_EACH_ID = False # imputation within each icustay_id\n",
    "IMPUTE_COLUMN = False # imputation based on whole column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data preprocessing\n",
    "from datetime import datetime\n",
    "print(datetime.now())\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import preprocessing\n",
    "from sklearn.model_selection import train_test_split\n",
    "import collections\n",
    "from collections import defaultdict\n",
    "#import time\n",
    "from datetime import datetime\n",
    "\n",
    "#logistic regression model: train and evaluation and XGB \n",
    "from sklearn.linear_model import LogisticRegression\n",
    "import xgboost as xgb\n",
    "from  sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "\n",
    "from sklearn.calibration import CalibratedClassifierCV, calibration_curve\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import roc_curve, precision_recall_curve, f1_score, roc_auc_score, auc, accuracy_score\n",
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.impute import SimpleImputer\n",
    "import matplotlib.lines as mlines\n",
    "import matplotlib.transforms as mtransforms\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the full files pathes are here\n",
    "# time dependent\n",
    "DATA_PATH_stages=\"../data/kdigo_stages_measured.csv\" \n",
    "DATA_PATH_labs = \"../data/labs-kdigo_stages_measured.csv\" \n",
    "DATA_PATH_vitals = \"../data/vitals-kdigo_stages_measured.csv\" \n",
    "DATA_PATH_vents = \"../data/vents-vasopressor-sedatives-kdigo_stages_measured.csv\"\n",
    "#no time dependent\n",
    "DATA_PATH_detail=\"../data/icustay_detail-kdigo_stages_measured.csv\" #age constraint\n",
    "#DATA_PATH_icd = \"../data/diagnoses_icd_aki_measured.csv\" #AL was \"...measured 2.csv\"\n",
    "\n",
    "SEPARATOR=\";\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameter as constant \n",
    "\n",
    "#which classifier to use, only run one classifier at one time \n",
    "ALL_STAGES = False\n",
    "CLASS1 = True   #AnyAKI\n",
    "#CLASS2 = True    #ModerateSevereAKI\n",
    "#CLASS3 = False    #SevereAKI\n",
    "\n",
    "MAX_DAYS = 35\n",
    "\n",
    "NORMALIZATION = 'min-max' \n",
    "\n",
    "CAPPING_THRESHOLD_UPPER = 0.99\n",
    "CAPPING_THRESHOLD_LOWER = 0.01\n",
    "\n",
    "#Age constraints: adults\n",
    "ADULTS_MIN_AGE = 18\n",
    "ADULTS_MAX_AGE = -1\n",
    "\n",
    "\n",
    "SPLIT_SIZE = 0.2\n",
    "\n",
    "#Two options to deal with time series data\n",
    "FIRST_TURN_POS = True #True #False # the first charttime that turn pos as bechmark for Target\n",
    "LAST_CHARTTIME = False #the last charttime as bechmark for Target\n",
    "\n",
    "# resampling or not\n",
    "TIME_SAMPLING = True \n",
    "SAMPLING_INTERVAL = '6H'\n",
    "RESAMPLE_LIMIT = 16 # 4 days*6h interval\n",
    "MOST_COMMON = False #resampling with most common\n",
    "# if MOST_COMMON is not applied,sampling with different strategies per kind of variable, \n",
    "# numeric variables use mean value, categorical variables use max value\n",
    "\n",
    "IMPUTE_METHOD = 'most_frequent' \n",
    "FILL_VALUE = 0 #fill missing value and ragged part of 3d array\n",
    "\n",
    "# How much time the prediction should occur (hours)\n",
    "HOURS_AHEAD = 48\n",
    "\n",
    "# precentage of dataset to be used for testing model\n",
    "TEST_SIZE = 0.1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set changable info corresponding to each classifier as variables\n",
    "\n",
    "min_set =  [\"icustay_id\", \"charttime\", \"creat\", \"uo_rt_6hr\", \"uo_rt_12hr\", \"uo_rt_24hr\", \"aki_stage\"]\n",
    "\n",
    "max_set = ['icustay_id', 'charttime', 'aki_stage', 'hadm_id','aniongap_avg', 'bicarbonate_avg', \n",
    "           'bun_avg','chloride_avg', 'creat', 'diasbp_mean', 'glucose_avg', 'heartrate_mean',\n",
    "           'hematocrit_avg', 'hemoglobin_avg', 'potassium_avg', 'resprate_mean','sodium_avg', 'spo2_mean', 'sysbp_mean', \n",
    "           'uo_rt_12hr', 'uo_rt_24hr','uo_rt_6hr', 'wbc_avg', 'sedative', 'vasopressor', 'vent', 'age', 'F','M', \n",
    "           'asian', 'black', 'hispanic', 'native', 'other', 'unknown','white', 'ELECTIVE', 'EMERGENCY', 'URGENT']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some functions used later\n",
    "\n",
    "\n",
    "def cap_data(df):\n",
    "    print(\"Capping between the {} and {} quantile\".format(CAPPING_THRESHOLD_LOWER, CAPPING_THRESHOLD_UPPER))\n",
    "    cap_mask = df.columns.difference(['icustay_id', 'charttime', 'aki_stage'])\n",
    "    df[cap_mask] = df[cap_mask].clip(df[cap_mask].quantile(CAPPING_THRESHOLD_LOWER),\n",
    "                                     df[cap_mask].quantile(CAPPING_THRESHOLD_UPPER),\n",
    "                                     axis=1)\n",
    "\n",
    "    return df\n",
    "\n",
    "# impute missing value in resampleing data with most common based on each id\n",
    "def fast_mode(df, key_cols, value_col):\n",
    "    \"\"\" Calculate a column mode, by group, ignoring null values. \n",
    "    \n",
    "    key_cols : list of str - Columns to groupby for calculation of mode.\n",
    "    value_col : str - Column for which to calculate the mode. \n",
    "\n",
    "    Return\n",
    "    pandas.DataFrame\n",
    "        One row for the mode of value_col per key_cols group. If ties, returns the one which is sorted first. \"\"\"\n",
    "    return (df.groupby(key_cols + [value_col]).size() \n",
    "              .to_frame('counts').reset_index() \n",
    "              .sort_values('counts', ascending=False) \n",
    "              .drop_duplicates(subset=key_cols)).drop('counts',axis=1)\n",
    "\n",
    "def normalise_data(df, norm_mask):\n",
    "    print(\"Normalizing in [0,1] with {} normalization\".format(NORMALIZATION))\n",
    "    \n",
    "    df[norm_mask] = (df[norm_mask] - df[norm_mask].min()) / (df[norm_mask].max() - df[norm_mask].min())\n",
    "    \n",
    "    return df\n",
    "\n",
    "\n",
    "#get max shape of 3d array\n",
    "def get_dimensions(array, level=0):   \n",
    "    yield level, len(array)\n",
    "    try:\n",
    "        for row in array:\n",
    "            yield from get_dimensions(row, level + 1)\n",
    "    except TypeError: #not an iterable\n",
    "        pass\n",
    "\n",
    "def get_max_shape(array):\n",
    "    dimensions = defaultdict(int)\n",
    "    for level, length in get_dimensions(array):\n",
    "        dimensions[level] = max(dimensions[level], length)\n",
    "    return [value for _, value in sorted(dimensions.items())]\n",
    "\n",
    "#pad the ragged 3d array to rectangular shape based on max size\n",
    "def iterate_nested_array(array, index=()):\n",
    "    try:\n",
    "        for idx, row in enumerate(array):\n",
    "            yield from iterate_nested_array(row, (*index, idx)) \n",
    "    except TypeError: # final level            \n",
    "        yield (*index, slice(len(array))), array # think of the types\n",
    "\n",
    "def pad(array, fill_value):\n",
    "    dimensions = get_max_shape(array)\n",
    "    result = np.full(dimensions, fill_value, dtype = np.float64)  \n",
    "    for index, value in iterate_nested_array(array):\n",
    "        result[index] = value \n",
    "    return result\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"read csv files\")\n",
    "#reading csv files\n",
    "X = pd.read_csv(DATA_PATH_stages, sep= SEPARATOR)\n",
    "X.drop([\"aki_stage_creat\", \"aki_stage_uo\"], axis = 1, inplace = True)\n",
    "#remove totally empty rows \n",
    "X = X.dropna(how = 'all', subset = ['creat','uo_rt_6hr','uo_rt_12hr','uo_rt_24hr','aki_stage'])\n",
    "print(\"convert charttime to timestamp\")\n",
    "X['charttime'] = pd.to_datetime(X['charttime'])\n",
    "#merge rows if they have exact timestamp within same icustay_id AL: t substitutes missing values with zero\n",
    "#X = X.groupby(['icustay_id', 'charttime']).sum().reset_index(['icustay_id', 'charttime'])\n",
    "\n",
    "dataset_detail = pd.read_csv(DATA_PATH_detail, sep= SEPARATOR)  #age constraint\n",
    "dataset_detail.drop(['dod', 'admittime','dischtime', 'los_hospital','ethnicity','hospital_expire_flag', 'hospstay_seq',\n",
    "       'first_hosp_stay', 'intime', 'outtime', 'los_icu', 'icustay_seq','first_icu_stay'], axis = 1, inplace = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_labs = pd.read_csv(DATA_PATH_labs, sep= SEPARATOR) # 'bands lactate platelet ptt inr pt\n",
    "dataset_labs.drop(['albumin_min', 'albumin_max','bilirubin_min', 'bilirubin_max','bands_min', 'bands_max',\n",
    "                   'lactate_min', 'lactate_max','platelet_min', 'platelet_max','ptt_min', 'ptt_max', \n",
    "                   'inr_min', 'inr_max', 'pt_min', 'pt_max'], axis = 1, inplace = True)\n",
    "dataset_labs = dataset_labs.dropna(subset=['charttime'])\n",
    "dataset_labs = dataset_labs.dropna(subset=dataset_labs.columns[4:], how='all')\n",
    "dataset_labs['charttime'] = pd.to_datetime(dataset_labs['charttime'])\n",
    "dataset_labs = dataset_labs.sort_values(by=['icustay_id', 'charttime'])\n",
    "\n",
    "if  MAX_FEATURE_SET:\n",
    "    dataset_vitals = pd.read_csv(DATA_PATH_vitals, sep= SEPARATOR)  \n",
    "    dataset_vents = pd.read_csv(DATA_PATH_vents , sep= SEPARATOR)\n",
    "    #dataset_icd = pd.read_csv(DATA_PATH_icd, sep= SEPARATOR)\n",
    "\n",
    "    dataset_vitals.drop([\"heartrate_min\", \"heartrate_max\",\"sysbp_min\", \"sysbp_max\",\"diasbp_min\", \"diasbp_max\",\n",
    "                        'meanbp_min','meanbp_max', 'meanbp_mean','tempc_min', 'tempc_max', 'tempc_mean',\n",
    "                        \"resprate_min\", \"resprate_max\", \"spo2_min\", \"spo2_max\", \"glucose_min\", \"glucose_max\"], axis = 1, inplace = True)\n",
    "          \n",
    "    print(\"convert charttime to timestamp\")\n",
    "    dataset_vitals['charttime'] = pd.to_datetime(dataset_vitals['charttime'])\n",
    "    dataset_vents['charttime'] = pd.to_datetime(dataset_vents['charttime'])\n",
    "    \n",
    "    dataset_vitals = dataset_vitals.sort_values(by=['icustay_id', 'charttime'])\n",
    "    dataset_vents = dataset_vents.sort_values(by=['icustay_id', 'charttime'])\n",
    "    \n",
    "    # AL drop those where all columns are nan\n",
    "    dataset_vitals = dataset_vitals.dropna(subset=dataset_vitals.columns[4:], how='all')   \n",
    "     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Labs file: instead of min and max their avg\n",
    "counter = 0\n",
    "col1 = 4\n",
    "col2 = 5\n",
    "null_l = [] # no null values in those that are different\n",
    "changed = 0 # 4316 records changed to avg\n",
    "\n",
    "while counter < 11:\n",
    "    row = 0\n",
    "# find where min and max are different and save their row indices \n",
    "    while row < len(dataset_labs):\n",
    "        a = dataset_labs.iloc[row,col1]\n",
    "        b = dataset_labs.iloc[row,col2]\n",
    "        if a==b or (np.isnan(a) and np.isnan(b)):\n",
    "            pass\n",
    "        elif a!=b:\n",
    "            changed +=1\n",
    "            avg = (a+b)/2\n",
    "            dataset_labs.iloc[row,col1] = avg\n",
    "            if (np.isnan(a) and ~np.isnan(b)) or (np.isnan(b) and ~np.isnan(a)):\n",
    "                null_l.append(row)\n",
    "        else:\n",
    "            print(a)\n",
    "            print(b)\n",
    "        row +=1       \n",
    "    # delete the redundant column max, update counters\n",
    "    dataset_labs.drop(dataset_labs.columns[col2], axis=1, inplace = True)\n",
    "    counter = counter+1\n",
    "    col1 = col1+1\n",
    "    col2 = col2+1\n",
    "\n",
    "dataset_labs.columns = ['subject_id','hadm_id', 'icustay_id', 'charttime', 'aniongap_avg', 'bicarbonate_avg', \n",
    "                        'creatinine_avg', 'chloride_avg', 'glucose_avg', 'hematocrit_avg','hemoglobin_avg',\n",
    "                        'potassium_avg', 'sodium_avg', 'bun_avg', 'wbc_avg']\n",
    "if len(null_l)>0:\n",
    "    print(\"null values encountered\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Merge creatinine and glucose.\")\n",
    "# merge creatinine from labs and set with labels\n",
    "creat_l = dataset_labs[['icustay_id','charttime','creatinine_avg']].copy()\n",
    "creat_l = creat_l.dropna(subset=['creatinine_avg'])\n",
    "creat = X[['icustay_id','charttime', 'creat']].copy()\n",
    "creat = creat.dropna(subset=['creat'])\n",
    "creat_l = creat_l.rename(columns={\"creatinine_avg\": \"creat\"})\n",
    "creat = creat.append(creat_l, ignore_index=True)\n",
    "creat.drop_duplicates(inplace = True)\n",
    "#delete old columns\n",
    "dataset_labs.drop([\"creatinine_avg\"], axis = 1, inplace = True)\n",
    "dataset_labs = dataset_labs.dropna(subset=dataset_labs.columns[4:], how='all')\n",
    "X.drop([\"creat\"], axis = 1, inplace = True)\n",
    "#merge new column\n",
    "X = pd.merge(X, creat, on = [\"icustay_id\", \"charttime\"], sort = True, how= \"outer\", copy = False)\n",
    "\n",
    "if MAX_FEATURE_SET:\n",
    "    # merge glucose from vitals and labs\n",
    "    glucose_v = dataset_vitals[['subject_id','hadm_id','icustay_id','charttime', 'glucose_mean']].copy()\n",
    "    glucose_v = glucose_v.dropna(subset=['glucose_mean'])\n",
    "    glucose = dataset_labs[['subject_id','hadm_id','icustay_id','charttime', 'glucose_avg']].copy()\n",
    "    glucose = glucose.dropna(subset=['glucose_avg'])\n",
    "    glucose_v = glucose_v.rename(columns={\"glucose_mean\": \"glucose_avg\"})\n",
    "    glucose = glucose.append(glucose_v, ignore_index=True)\n",
    "    glucose.drop_duplicates(inplace = True)\n",
    "    #delete old columns\n",
    "    dataset_labs.drop([\"glucose_avg\"], axis = 1, inplace = True)\n",
    "    dataset_vitals.drop([\"glucose_mean\"], axis = 1, inplace = True)\n",
    "    dataset_vitals = dataset_vitals.dropna(subset=dataset_vitals.columns[4:], how='all')\n",
    "    #merge new column\n",
    "    dataset_labs = pd.merge(dataset_labs, glucose, on = ['subject_id','hadm_id','icustay_id','charttime',], sort = True, how= \"outer\", copy = False)\n",
    "    \n",
    "dataset_labs = dataset_labs.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)\n",
    "X = X.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Merging labs, vitals and vents files\")\n",
    "#merge files with time-dependent data, based on icustay_id and charttime\n",
    "if MAX_FEATURE_SET:\n",
    "    X = pd.merge(X, dataset_labs, on = [\"icustay_id\", \"charttime\"], how= \"outer\", copy = False)\n",
    "    X = pd.merge(X, dataset_vitals, on = [\"icustay_id\", \"charttime\",\"subject_id\", \"hadm_id\"], how= \"outer\", copy = False)\n",
    "    X = pd.merge(X, dataset_vents, on = [\"icustay_id\", \"charttime\"], how= \"outer\", copy = False) \n",
    "    X.drop([\"subject_id\"], axis = 1, inplace = True) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"start preprocessing time dependent data\") \n",
    "print(\"Removing patients under the min age\")\n",
    "dataset_detail = dataset_detail.loc[dataset_detail['age'] >= ADULTS_MIN_AGE]\n",
    "adults_icustay_id_list = dataset_detail['icustay_id'].unique()\n",
    "X = X[X.icustay_id.isin(adults_icustay_id_list)].sort_values(by=['icustay_id'], ignore_index = True)\n",
    "X = X.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)\n",
    "adults_icustay_id_list = np.sort(adults_icustay_id_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"drop icustay_id with time span less than 48hrs\")\n",
    "def more_than_HOURS_ahead(adults_icustay_id_list, X):\n",
    "    drop_list = []\n",
    "    los_list = [] # calculating LOS ICU based on charttime\n",
    "    long_stays_id = [] # LOS longer than MAX DAYS days\n",
    "    last_charttime_list = []\n",
    "    seq_length = X.groupby(['icustay_id'],as_index=False).size()\n",
    "    id_count = 0\n",
    "    first_row_index = 0\n",
    "\n",
    "    while id_count < len(adults_icustay_id_list):\n",
    "        icustay_id = adults_icustay_id_list[id_count]\n",
    "        last_row_index = first_row_index + seq_length.iloc[id_count,1]-1\n",
    "        first_time = X.iat[first_row_index, X.columns.get_loc('charttime')]\n",
    "        last_time = X.iat[last_row_index, X.columns.get_loc('charttime')]\n",
    "        los = round(float((last_time - first_time).total_seconds()/60/60/24),4) # in days\n",
    "        if los < HOURS_AHEAD/24:\n",
    "            drop_list.append(icustay_id)\n",
    "        else:\n",
    "            los_list.append(los)\n",
    "            if los > MAX_DAYS:\n",
    "                long_stays_id.append(icustay_id)\n",
    "                last_charttime_list.append(last_time)\n",
    "        # udpate for the next icustay_id\n",
    "        first_row_index = last_row_index+1\n",
    "        id_count +=1\n",
    "    if len(long_stays_id) != len(last_charttime_list):\n",
    "        print('ERROR')\n",
    "    print(\"%d long stays\" % len(long_stays_id))\n",
    "    # drop all the rows with the saved icustay_id\n",
    "    print(\"there are %d id-s shorter than 48 hours\" % len(drop_list))\n",
    "    X = X[~X.icustay_id.isin(drop_list)]\n",
    "    id_list = X['icustay_id'].unique()\n",
    "    X = X.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)\n",
    "    \n",
    "    return id_list, X, long_stays_id,last_charttime_list\n",
    "\n",
    "id_list, X, long_stays_id,last_charttime_list  = more_than_HOURS_ahead(adults_icustay_id_list, X)\n",
    "\n",
    "\n",
    "long = pd.DataFrame()\n",
    "long['icustay_id']  = long_stays_id\n",
    "long['last_time']  = last_charttime_list\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# deleting rows that are not within MAX_DAYS (35) period\n",
    "i = 0 # long df index\n",
    "drop_long_time = []\n",
    "    \n",
    "while i < len(long_stays_id):\n",
    "    j = 0\n",
    "    all_rows = X.index[X['icustay_id'] == long.loc[i,'icustay_id']].tolist()\n",
    "    while j < len(all_rows):\n",
    "        time = X.iat[all_rows[j], X.columns.get_loc('charttime')]\n",
    "        # if keep last MAX_DAYS \n",
    "        if (long.loc[i,'last_time'] - time).total_seconds() > MAX_DAYS*24*60*60:\n",
    "            drop_long_time.append(all_rows[j])\n",
    "            j +=1\n",
    "        else:\n",
    "            break\n",
    "    i +=1       \n",
    "X.drop(X.index[drop_long_time], inplace=True) \n",
    "\n",
    "# checking for 48h min length again\n",
    "id_list, X, long_stays_id,last_charttime_list  = more_than_HOURS_ahead(id_list, X)\n",
    "dataset_detail = dataset_detail[dataset_detail.icustay_id.isin(id_list)].sort_values(by=['icustay_id'], ignore_index = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For testing purpose, use small amount of data first\n",
    "if TESTING:\n",
    "    rest, id_list = train_test_split(id_list, test_size= 0.01, random_state=42)\n",
    "    X = X[X.icustay_id.isin(id_list)].sort_values(by=['icustay_id'])\n",
    "    dataset_detail = dataset_detail[dataset_detail.icustay_id.isin(id_list)].sort_values(by=['icustay_id'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resample  and impute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (TIME_SAMPLING and MOST_COMMON):\n",
    "    print(\"resampling: MOST_COMMON\")\n",
    "    # Resample the data using assigned interval,mode() for most common\n",
    "    X = X.set_index('charttime').groupby('icustay_id').resample(SAMPLING_INTERVAL).mode().reset_index()\n",
    "elif TIME_SAMPLING:\n",
    "    print(\"resampling: MEAN & ZERO\")\n",
    "    # Sampling with different strategies per kind of variable\n",
    "    label = ['aki_stage']\n",
    "    skip = ['icustay_id', 'charttime', 'aki_stage']\n",
    "    if MAX_FEATURE_SET:\n",
    "        discrete_feat = ['sedative', 'vasopressor', 'vent', 'hadm_id']\n",
    "        skip.extend(discrete_feat)    \n",
    "    # all features that are not in skip are numeric\n",
    "    numeric_feat = list(X.columns.difference(skip))\n",
    "    \n",
    "    # Applying aggregation to features depending on their type\n",
    "    X = X.set_index('charttime').groupby('icustay_id').resample(SAMPLING_INTERVAL)\n",
    "    if MAX_FEATURE_SET:\n",
    "        X_discrete = X[discrete_feat].max().fillna(FILL_VALUE).astype(np.int64)\n",
    "    X_numeric = X[numeric_feat].mean() \n",
    "    X_label = X['aki_stage'].max()\n",
    "    print(\"Merging sampled features\")\n",
    "    try:\n",
    "        X = pd.concat([X_numeric, X_discrete,X_label], axis=1).reset_index()\n",
    "    except:\n",
    "        X = pd.concat([X_numeric,X_label], axis=1).reset_index()\n",
    "print(X.shape)\n",
    "#Label forward fill\n",
    "X['aki_stage'] = X['aki_stage'].ffill(limit=RESAMPLE_LIMIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Imputation.\")\n",
    "# do imputation of label with zero if there are still missing values\n",
    "X['aki_stage'] = X['aki_stage'].fillna(0)\n",
    "# using most common within each icustay_id\n",
    "if IMPUTE_EACH_ID:\n",
    "    # set a new variable so won't change the orginial X\n",
    "    column_name = list(X.columns)\n",
    "    column_name.remove(column_name[0]) \n",
    "    for feature in column_name:\n",
    "        X.loc[X[feature].isnull(), feature] = X.icustay_id.map(fast_mode(X, ['icustay_id'], feature).set_index('icustay_id')[feature])       \n",
    "\n",
    "# imputation based on whole column\n",
    "if IMPUTE_COLUMN:\n",
    "    imp = SimpleImputer(missing_values=np.nan, strategy= IMPUTE_METHOD)\n",
    "    cols = list(X.columns)\n",
    "    cols = cols[2:23]\n",
    "    X[cols]=imp.fit_transform(X[cols])  \n",
    "\n",
    "# If no imputation method selected or only impute each id, for the remaining nan impute direclty with FILL_VALUE\n",
    "X = X.fillna(FILL_VALUE) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# more comfortable to review in this order\n",
    "try:\n",
    "    cols = ['icustay_id', 'charttime','aki_stage','hadm_id','aniongap_avg','bicarbonate_avg', 'bun_avg','chloride_avg',\n",
    "            'creat','diasbp_mean', 'glucose_avg', 'heartrate_mean', 'hematocrit_avg','hemoglobin_avg', \n",
    "            'potassium_avg', 'resprate_mean', 'sodium_avg','spo2_mean', 'sysbp_mean', 'uo_rt_12hr', \n",
    "            'uo_rt_24hr', 'uo_rt_6hr','wbc_avg', 'sedative', 'vasopressor', 'vent' ]\n",
    "    X = X[cols]\n",
    "    print(\"success\")\n",
    "except:\n",
    "    try:\n",
    "        cols = ['icustay_id', 'charttime','aki_stage','creat','uo_rt_12hr', 'uo_rt_24hr', 'uo_rt_6hr']\n",
    "        X = X[cols]\n",
    "    except:\n",
    "        print(\"error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"binarise labels\")\n",
    "if ALL_STAGES:\n",
    "    pass\n",
    "elif CLASS1:\n",
    "    X.loc[X['aki_stage'] > 1, 'aki_stage'] = 1\n",
    "elif CLASS2:\n",
    "    X.loc[X['aki_stage'] < 2, 'aki_stage'] = 0\n",
    "    X.loc[X['aki_stage'] > 1, 'aki_stage'] = 1\n",
    "elif CLASS3:\n",
    "    X.loc[X['aki_stage'] < 3, 'aki_stage'] = 0\n",
    "    X.loc[X['aki_stage'] > 2, 'aki_stage'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"Shifting the labels 48 h\") # by 8 position : 6h sampling*8=48h and ffil 8 newly empty ones\n",
    "# group by\n",
    "X['aki_stage'] = X.groupby('icustay_id')['aki_stage'].shift(-(HOURS_AHEAD // int(SAMPLING_INTERVAL[:-1])))\n",
    "X = X.dropna(subset=['aki_stage'])\n",
    "X['icustay_id'].nunique()\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if LAST_CHARTTIME:\n",
    "    print(\"Here we choose the last charttime as bechmark for Target\")\n",
    "    print(\"find last charttime for each icustay_id\")\n",
    "    #extract the aki_stage value of row of \"last\" charttime of each icustay_id as the label/target set\n",
    "    #find last charttime for each icustay_id\n",
    "\n",
    "    target_list = [] #to store target(y)\n",
    "    last_charttime_list= [] #a list to store last charttime for every icustay_id\n",
    "    last_charttime_index_list = [] #a list to store index when last charttime for every icustay_id\n",
    "\n",
    "    dataset_size = X.shape[0]\n",
    "    k=0 #index\n",
    "    h=0 #the times that id matches\n",
    "\n",
    "    for extract_id in id_list:\n",
    "        while k < X:\n",
    "            if X.iat[k, X.columns.get_loc('icustay_id')] == extract_id:\n",
    "                if h==0 and k != dataset_size-1:\n",
    "                    k = k+1\n",
    "                    h = h+1\n",
    "                elif h !=0:\n",
    "                    if k == dataset_size-1:\n",
    "                        target_list.append(X.iat[k, X.columns.get_loc('aki_stage')])\n",
    "                        last_charttime_list.append(X.iat[k, X.columns.get_loc('charttime')])\n",
    "                        last_charttime_index_list.append(k)\n",
    "                    k = k+1\n",
    "            elif X.iat[k, X.columns.get_loc('icustay_id')] != extract_id:\n",
    "                target_list.append(X.iat[k-1, X.columns.get_loc('aki_stage')])\n",
    "                last_charttime_list.append(X.iat[k-1, X.columns.get_loc('charttime')])\n",
    "                last_charttime_index_list.append(k-1)\n",
    "                h=0\n",
    "                break\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if FIRST_TURN_POS:\n",
    "    # create one label per icustay_id - first turn positive approach\n",
    "\n",
    "    X = X.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)\n",
    "    id_list.sort()\n",
    "    last_charttime_list= []\n",
    "\n",
    "    index_list = []\n",
    "    label_list = []\n",
    "\n",
    "    first_row_index = 0\n",
    "    id_count = 0\n",
    "    seq_length = X.groupby(['icustay_id'],as_index=False).size()\n",
    "\n",
    "    for ID in id_list:\n",
    "        last_row_index = first_row_index + seq_length.iloc[id_count,1]-1\n",
    "        a = X.loc[X['icustay_id']==ID].aki_stage\n",
    "        if 1 not in a.values:\n",
    "            label_list.append(0)\n",
    "            last_charttime_list.append(X.iat[last_row_index, X.columns.get_loc('charttime')]) \n",
    "            index_list.append(last_row_index)\n",
    "        elif 1 in a.values:\n",
    "            label_list.append(1)\n",
    "            row = first_row_index\n",
    "            while row != last_row_index+1:\n",
    "                if X.iat[row, X.columns.get_loc('aki_stage')]==0:\n",
    "                    row +=1\n",
    "                elif X.iat[row, X.columns.get_loc('aki_stage')]==1:\n",
    "                    last_charttime_list.append(X.iat[row, X.columns.get_loc('charttime')])\n",
    "                    index_list.append(row)\n",
    "                    break\n",
    "        first_row_index = last_row_index+1\n",
    "        id_count +=1 \n",
    "    X = X.drop(['aki_stage'], axis=1)\n",
    "    Thresholds = pd.DataFrame({'icustay_id':id_list, 'charttime': last_charttime_list, 'final_label':label_list})\n",
    "    X = (Thresholds.merge(X, on='icustay_id', how='left',suffixes=('_x','')).query(\"charttime_x >= charttime\").reindex(columns=X.columns))        \n",
    "    print(X.shape)\n",
    "    X['icustay_id'].nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_register =  pd.DataFrame()\n",
    "label_register['icustay_id'] = id_list\n",
    "label_register['label'] = label_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Add categorical features (details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"start preprocessing not time dependent data\")\n",
    "if MAX_FEATURE_SET:\n",
    "    #extract datasets based on id_list\n",
    "    dataset_detail = dataset_detail.loc[dataset_detail['icustay_id'].isin(id_list)]\n",
    "    #sort by ascending order\n",
    "    dataset_detail = dataset_detail.sort_values(by=['icustay_id'])\n",
    "    subject_id = dataset_detail[\"subject_id\"].unique()\n",
    "    #transfrom categorical data to binary form\n",
    "    dataset_detail = dataset_detail.join(pd.get_dummies(dataset_detail.pop('gender')))\n",
    "    dataset_detail = dataset_detail.join(pd.get_dummies(dataset_detail.pop(\"ethnicity_grouped\")))\n",
    "    dataset_detail = dataset_detail.join(pd.get_dummies(dataset_detail.pop('admission_type')))\n",
    "    dataset_detail = dataset_detail.drop(['subject_id', 'hadm_id'], axis=1)\n",
    "    X =  pd.merge(X, dataset_detail, on = [\"icustay_id\"], how= \"left\", copy = False) \n",
    "    numeric_feat.append('age')\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cap features between 0.01 / 0.99 quantile and normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = cap_data(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = normalise_data(X, numeric_feat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set target as y\n",
    "'''\n",
    "if LAST_CHARTTIME:\n",
    "    y = pd.DataFrame(target_list)\n",
    "if FIRST_TURN_POS:\n",
    "    y = Thresholds['final_label']\n",
    "print(len(y))\n",
    "y = np.array(y)\n",
    "'''\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"divide dataset into train, test and validation sets\")\n",
    "id_train, id_test_val = train_test_split(id_list, test_size = SPLIT_SIZE, random_state = 42) # train set is 80%)\n",
    "# to ensure the same split as LSTM\n",
    "id_valid, id_test = train_test_split(id_test_val, test_size = 0.5, random_state = 42) # test 10% valid 10%\n",
    "print(\" test is %d\" %len(id_test))\n",
    "id_train = np.concatenate((id_train, id_valid), axis=None)\n",
    "print(\"train is %d\" % len(id_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X[X.icustay_id.isin(id_train)].sort_values(by=['icustay_id'])\n",
    "X_test = X[X.icustay_id.isin(id_test)].sort_values(by=['icustay_id'], ignore_index = True) \n",
    "X_test = X_test.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)\n",
    "X_train = X_train.sort_values(by=['icustay_id', 'charttime'], ignore_index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X_train.drop(['hadm_id', 'charttime'], axis=1)\n",
    "X_test = X_test.drop(['hadm_id', 'charttime'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test = []\n",
    "id_ordered = X_test['icustay_id'].unique()\n",
    "for i in id_ordered:\n",
    "    a = label_register.loc[label_register['icustay_id']==i].label\n",
    "    y_test.append(a.values[0])\n",
    "    \n",
    "y_train = []\n",
    "id_ordered = X_train['icustay_id'].unique()\n",
    "for i in id_ordered:\n",
    "    a = label_register.loc[label_register['icustay_id']==i].label\n",
    "    y_train.append(a.values[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"reshape 2D dataframe to 3D Array, group by icustay_id\")\n",
    "X_train = np.array(list(X_train.groupby(['icustay_id'],as_index=False).apply(pd.DataFrame.to_numpy)))\n",
    "print(X_train.shape)\n",
    "#print(\"reshape 2D dataframe to 3D Array, group by icustay_id\")\n",
    "X_test = np.array(list(X_test.groupby(['icustay_id'],as_index=False).apply(pd.DataFrame.to_numpy)))\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"fill ragged part(missing value) of 3d array with FILL_VALUE\")\n",
    "X_train = pad(X_train, FILL_VALUE)\n",
    "X_train = X_train[:,:,1:] #delete icustay id\n",
    "print(X_train.shape)\n",
    "X_test = pad(X_test, FILL_VALUE)\n",
    "X_test = X_test[:,:,1:] #delete icustay id\n",
    "print(X_test.shape)\n",
    "features = X_train.shape[2]\n",
    "timestamps = X_train.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"flatten 3D array back to 2D array\")\n",
    "nsamples, nx, ny = X_train.shape\n",
    "X_train = X_train.reshape((nsamples,nx*ny))\n",
    "print(X_train.shape)\n",
    "nsamples, nx, ny = X_test.shape\n",
    "X_test = X_test.reshape((nsamples,nx*ny))\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.amax(X_train) #test\n",
    "np.amax(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def performance (y_test, pred_probabilities):\n",
    "    # performance\n",
    "    fpr, tpr, thresholds = roc_curve(y_test, pred_probabilities)\n",
    "    # compute roc auc\n",
    "    roc_auc = roc_auc_score(y_test, pred_probabilities, average = 'micro')\n",
    "    # compute Precision_Recall curves\n",
    "    precision, recall, _ = precision_recall_curve(y_test, pred_probabilities)\n",
    "    # compute PR_AUC\n",
    "    pr_auc = metrics.auc(recall, precision)   \n",
    "    # I add confusion matrix\n",
    "    optimal_cut_off = round(thresholds[np.argmax(tpr - fpr)],4)\n",
    "    a = np.where(pred_probabilities > optimal_cut_off, 1, 0)\n",
    "    brier = round(metrics.brier_score_loss(y_test, pred_probabilities, sample_weight=None, pos_label=None),3)\n",
    "    predictions = np.where(pred_probabilities > optimal_cut_off, 1, 0)  \n",
    "    \n",
    "    print (\"Area Under ROC Curve: %0.2f\" % roc_auc  )\n",
    "    #print (\"Area Under PR Curve(AP): %0.2f\" % pr_auc  ) \n",
    "    print(\"Brier score : {:.3f}\".format(brier))\n",
    "    #print('Accuracy for Classifier : {:.2f}'.format(accuracy_score(y_test, predictions)))\n",
    "    #print('Cut off: ' + str(optimal_cut_off))\n",
    "    matrix = metrics.confusion_matrix(y_test, a, labels=None, normalize=None)\n",
    "    #print(str(matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bin_total(y_true, y_prob, n_bins):\n",
    "    bins = np.linspace(0., 1. + 1e-8, n_bins + 1)\n",
    "\n",
    "    # In sklearn.calibration.calibration_curve, the last value in the array is always 0.\n",
    "    binids = np.digitize(y_prob, bins) - 1\n",
    "\n",
    "    return np.bincount(binids, minlength=len(bins))\n",
    "\n",
    "def missing_bin(bin_array):\n",
    "    midpoint = \" \"    \n",
    "    if bin_array[0]==0:\n",
    "        midpoint = \"5%, \"\n",
    "    if bin_array[1]==0:\n",
    "        midpoint = midpoint + \"15%, \"\n",
    "    if bin_array[2]==0:\n",
    "        midpoint = midpoint + \"25%, \"\n",
    "    if bin_array[3]==0:\n",
    "        midpoint = midpoint + \"35%, \" \n",
    "    if bin_array[4]==0:\n",
    "        midpoint = midpoint + \"45%, \"\n",
    "    if bin_array[5]==0:\n",
    "        midpoint = midpoint + \"55%, \"\n",
    "    if bin_array[6]==0:\n",
    "        midpoint = midpoint + \"65%, \"\n",
    "    if bin_array[7]==0:\n",
    "        midpoint = midpoint + \"75%, \"\n",
    "    if bin_array[8]==0:\n",
    "        midpoint = midpoint + \"85%, \"\n",
    "    if bin_array[9]==0:\n",
    "        midpoint = midpoint + \"95%, \"\n",
    "    return \"The missing bins have midpoint values of \"+ str(midpoint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_graphs (y_test,pred_probabilities, classifier_name, plot_name, algorithm):\n",
    "    \n",
    "    # performance\n",
    "    fpr, tpr, thresholds = roc_curve(y_test, pred_probabilities)\n",
    "    # compute roc auc\n",
    "    roc_auc = roc_auc_score(y_test, pred_probabilities, average = 'micro')\n",
    "    # compute Precision_Recall curves\n",
    "    precision, recall, _ = precision_recall_curve(y_test, pred_probabilities)\n",
    "    # compute PR_AUC\n",
    "    pr_auc = metrics.auc(recall, precision)\n",
    "\n",
    "    # compute calibration curve\n",
    "    LR_y, LR_x = calibration_curve(y_test, pred_probabilities, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array = bin_total(y_test, pred_probabilities , n_bins=10)\n",
    "    print(missing_bin(bin_array))\n",
    "\n",
    "    print(\"plot curves and save in one png file\")\n",
    "    #save three plots in one png file\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(3, figsize=(7, 24))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    fig.suptitle('Evaluation of '+ plot_name)\n",
    "\n",
    "    fpr, tpr, thresholds = roc_curve(y_test, pred_probabilities)\n",
    "    \n",
    "    # plot roc curve\n",
    "    ax1.plot(fpr, tpr, 'C2', label=algorithm+\" \"+\"Classifier \" + str(classifier_name) + \", auc=\" +str(round(roc_auc,2)))\n",
    "    ax1.title.set_text('ROC AUC')\n",
    "    ax1.set(xlabel='False Positive Rate', ylabel='True Positive Rate')\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "\n",
    "    # plot PR curve\n",
    "    ax2.plot(recall, precision, 'C2', label=algorithm+\" \"+\"Classifier \" + str(classifier_name) + \", auc=\"+str(round(pr_auc,2)))\n",
    "    ax2.title.set_text('PR AUC')\n",
    "    ax2.set(xlabel='Recall', ylabel='Precision')\n",
    "    ax2.legend(loc=\"lower right\")\n",
    "\n",
    "    # plot calibration curve\n",
    "    ax3.plot(LR_x, LR_y, 'C2',marker='o', linewidth=1, label='LR')\n",
    "    line = mlines.Line2D([0, 1], [0, 1], color='black')\n",
    "    transform = ax3.transAxes\n",
    "    line.set_transform(transform)\n",
    "    ax3.add_line(line)\n",
    "    ax3.title.set_text('Calibration plot for '+str(plot_name))\n",
    "    ax3.set(xlabel= 'Predicted probability', ylabel= 'True probability in each bin')\n",
    "    ax3.legend(loc=\"lower right\")\n",
    "\n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distribution(pred_probabilities, dist_name):\n",
    "    #probabilities distributions graphs\n",
    "    true_1 = pd.DataFrame(pred_probabilities, columns=['Predicted probabilities'])\n",
    "    true_1['labels'] = y_test\n",
    "    true_0 = true_1.copy(deep = True) \n",
    "    print(true_1.shape)\n",
    "    indexNames = true_1[ true_1['labels'] == 0 ].index\n",
    "    true_1.drop(indexNames , inplace=True)\n",
    "    indexNames = true_0[ true_0['labels'] == 1 ].index\n",
    "    true_0.drop(indexNames , inplace=True)\n",
    "    true_1.drop(columns=['labels'], inplace = True)\n",
    "    true_0.drop(columns=['labels'], inplace = True)\n",
    "    \n",
    "    sns.distplot(true_1['Predicted probabilities'], hist = False, kde = True,\n",
    "                 kde_kws = {'shade': True, 'linewidth': 3, \"color\": \"g\"}, label = 'Class 1')\n",
    "    plt.ylabel('Density')\n",
    "    sns.distplot(true_0['Predicted probabilities'], hist = False, kde = True,\n",
    "                     kde_kws = {'shade': True, 'linewidth': 3}, label = 'Class 0')\n",
    "    plt.title('Density Plot'+ dist_name)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RUN MODELS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LR with best parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = ['Anion gap', 'Bicarbonate', 'Blood Urea Nitrogen', 'Chloride', 'Creatinine', 'Diastolic BP', 'Glucose', 'Heart rate', \n",
    "            'Hematocrit', 'Hemoglobin', 'Potassium', 'Respiratory rate', 'Sodium', 'Oxygen saturation', 'Systolic BP', 'Urine output 12h', 'Urine output 24h', 'Urine output 6h',\n",
    "            'White cell count', 'Sedative', 'Vasopressor', 'Ventilation', 'Age', 'Female gender', 'Male gender', 'Asian ethnicity', 'Black ethnicity', 'Hispanic ethnicity', 'Native american', \n",
    "            'Other ethnicity', 'Ethnicity unknown', 'White ethnicity', 'Elective admission', 'Emergency admission', 'Urgent admission']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(datetime.now())\n",
    "print(\"build LR\")\n",
    "model_LR1 = LogisticRegression(penalty = \"l1\", C = 0.1,fit_intercept = True, solver = 'saga', max_iter= 1000, \n",
    "                                  class_weight= 'balanced', random_state=42)\n",
    "\n",
    "print(\"train LR\")\n",
    "model_LR1.fit(X_train, y_train) \n",
    "# Make predictions\n",
    "LR1_pred = model_LR1.predict(X_test)\n",
    "# Predict probabilities on testing points\n",
    "LR1_probs = model_LR1.predict_proba(X_test)[::,1]\n",
    "print(datetime.now())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "performance(y_test,LR1_probs)\n",
    "\n",
    "'''\n",
    "# NO IMP \n",
    "L1\n",
    "Area Under ROC Curve: 0.90\n",
    "Brier score : 0.128\n",
    "\n",
    " \n",
    "# WITH IMP \n",
    "\n",
    "L1\n",
    "Area Under ROC Curve: 0.84\n",
    "Brier score : 0.161\n",
    "\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LogR Feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The estimated coefficients will all be around 1:\n",
    "LRcoefs = model_LR1.coef_\n",
    "#print(model_LR1.coef_)\n",
    "#LRcoefsStd = np.std(X, 0)*model_LR1.coef_\n",
    "#print(np.std(X, 0)*model_LR1.coef_)\n",
    "#LRcoefsStd = LRcoefsStd.flatten()\n",
    "\n",
    "LRcoefs = LRcoefs.flatten()\n",
    "LRcoefs.shape # 4655 = 133*35"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape to timestamps x features\n",
    "logR_feat = np.reshape(LRcoefs,(timestamps,-1)) # shape 133, 35\n",
    "# record the sign of coefficients\n",
    "pos_neg = logR_feat.mean(axis=0)\n",
    "signs = []\n",
    "for i in pos_neg:\n",
    "    if i<0:\n",
    "        signs.append(-1)\n",
    "    else:\n",
    "        signs.append(1)\n",
    "signs = np.array(signs)\n",
    "logR_feat = np.exp(logR_feat)\n",
    "logR_feat = logR_feat.mean(axis = 0)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0\n",
    "for i in signs:\n",
    "    logR_feat[a] = logR_feat[a]*i\n",
    "    a +=1\n",
    "logR_df =  pd.DataFrame(logR_feat, columns = ['Feature Importance'])\n",
    "logR_df['Features'] = variables\n",
    "logR_df = logR_df.sort_values(by = ['Feature Importance'], ascending = False, ignore_index = True)\n",
    "#logR_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.barplot(x='Feature Importance', y='Features', data=logR_df)\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 6)\n",
    "#plt.xlim(0.925, 1.1)\n",
    "plt.title('Logistic Regression feature Importance')\n",
    "plt.savefig('LogR_feature_importance_sign.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply signs\n",
    "a = 0\n",
    "for i in signs:\n",
    "    logR_feat[a] = logR_feat[a]*i\n",
    "    a +=1\n",
    "#print(logR_feat)    \n",
    "logR_df =  pd.DataFrame(logR_feat, columns = ['Feature Importance'])\n",
    "logR_df['Features'] = variables\n",
    "logR_df['abs'] = abs(logR_df['Feature Importance'])\n",
    "logR_df = logR_df.sort_values(by = ['abs'], ascending = False, ignore_index = True)\n",
    "logR_df_10 = logR_df.head(10)\n",
    "#logR_df_10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.barplot(x='Feature Importance', y='Features', data=logR_df_10, palette=\"mako\")\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 10)\n",
    "plt.xlim(0.99, 1.1)\n",
    "plt.title('Logistic Regression top 10 features by feature importance')\n",
    "plt.savefig('LogR_top10_feature_importance.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XGB model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counter=collections.Counter(y_train)\n",
    "SCALE_POS_WEIGHT = round(counter[0]/counter[1],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"build GBT\")\n",
    "#build gradient boost model\n",
    "xgb_model = xgb.XGBClassifier(\n",
    " learning_rate = 0.1,\n",
    " n_estimators= 70,\n",
    " max_depth= 7,\n",
    " min_child_weight= 3,\n",
    " gamma= 0.5,\n",
    " subsample= 0.8,\n",
    " colsample_bytree= 0.8,\n",
    " objective= 'binary:logistic',\n",
    " nthread= 4,\n",
    " scale_pos_weight= SCALE_POS_WEIGHT,\n",
    " seed= 42)\n",
    "\n",
    "print(\"train GBT\")\n",
    "# train xgb\n",
    "xgb_model.fit(X_train, y_train, \n",
    "                early_stopping_rounds= 40, \n",
    "                eval_metric= [\"auc\", \"error\", \"logloss\"], \n",
    "                eval_set= [(X_test, y_test)],\n",
    "                verbose= False)\n",
    "print(\"predictions after tuning\")\n",
    "#Make predictions \n",
    "GBT_pred = xgb_model.predict(X_test)\n",
    "GBT_probs = xgb_model.predict_proba(X_test)[::,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "performance(y_test,GBT_probs)\n",
    "'''\n",
    "No imp\n",
    "\n",
    "Area Under ROC Curve: 0.92\n",
    "Brier score : 0.108\n",
    " \n",
    " With imp\n",
    "Area Under ROC Curve: 0.91\n",
    "Brier score : 0.120\n",
    " '''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XGB feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape to timestamps x features\n",
    "xgb_feat = np.reshape(xgb_model.feature_importances_,(timestamps,-1))\n",
    "# sum the feature importance by feature\n",
    "xgb_feat = xgb_feat.sum(axis=0)\n",
    "xgb_feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_feat.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_df =  pd.DataFrame(xgb_feat, columns = ['Feature Importance'])\n",
    "xgb_df['Features'] = variables\n",
    "#xgb_df = xgb_df.set_index('Features')\n",
    "xgb_df = xgb_df.sort_values(by = ['Feature Importance'], ascending = False, ignore_index = True)\n",
    "#xgb_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.barplot(x='Feature Importance', y='Features', data=xgb_df)\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 6)\n",
    "plt.title('XGB feature Importance')\n",
    "plt.savefig('XGB_feature_importance.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_df_10 = xgb_df.head(10)\n",
    "ax = sns.barplot(x='Feature Importance', y='Features', data=xgb_df_10, palette = 'mako')\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 10)\n",
    "plt.title('XGB top 10 features by features importance')\n",
    "plt.savefig('XGB_top10_feature_importance.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Random Forests model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"build RF\")\n",
    "# building RF (default parameters)\n",
    "rf_model = RandomForestClassifier(n_estimators=700,criterion='gini',max_depth=40,min_samples_split=10,min_samples_leaf=2,\n",
    "                                  min_weight_fraction_leaf=0.0,max_features='auto',max_leaf_nodes=None,min_impurity_decrease=0.0,\n",
    "                                  min_impurity_split=None,bootstrap=False,oob_score=False,n_jobs=None,random_state=None,verbose=0,\n",
    "                                  warm_start=False,class_weight=None,ccp_alpha=0.0,max_samples=None)\n",
    "\n",
    "\n",
    "print(\"train RF\")\n",
    "# train RF\n",
    "rf_model.fit(X_train, y_train)\n",
    "\n",
    "print(\"predictions\")\n",
    "#Make predictions \n",
    "RF_pred = rf_model.predict(X_test)\n",
    "RF_probs = rf_model.predict_proba(X_test)[::,1] # probability of class 1\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "performance(y_test,RF_probs)\n",
    "'''\n",
    "No Imp\n",
    "\n",
    "Area Under ROC Curve: 0.93\n",
    "Brier score : 0.105\n",
    "\n",
    "With imp\n",
    "Area Under ROC Curve: 0.89\n",
    "Brier score : 0.133\n",
    " '''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RF feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape to timestamps x features\n",
    "rf_feat = np.reshape(rf_model.feature_importances_,(timestamps,-1))\n",
    "# sum the feature importance by feature\n",
    "rf_feat = rf_feat.sum(axis=0)\n",
    "#rf_feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_df =  pd.DataFrame(rf_feat, columns = ['Feature Importance'])\n",
    "rf_df['Features'] = variables\n",
    "#rf_df = rf_df.set_index('Features')\n",
    "rf_df = rf_df.sort_values(by = ['Feature Importance'], ascending = False, ignore_index = True)\n",
    "#rf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.barplot(x='Feature Importance', y='Features', data=rf_df)\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 6)\n",
    "plt.title('RF feature Importance')\n",
    "plt.savefig('RF_feature_importance.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_df_10 = rf_df.head(10)\n",
    "ax = sns.barplot(x='Feature Importance', y='Features', data=rf_df_10, palette = 'mako')\n",
    "ax.set_xlabel('Feature Importance', fontsize = 15)\n",
    "ax.set_ylabel(\"Features\",fontsize=15)\n",
    "ax.set_yticklabels(ax.get_ymajorticklabels(), fontsize = 10)\n",
    "plt.title('RF top 10 features by feature importance')\n",
    "plt.savefig('RF_top10_feature_importance.png', dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# calibration together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calibration_together (y_test,x1, x2,x3,x4, plot_name):\n",
    "    \n",
    "    # performance\n",
    "    #fpr, tpr, thresholds = roc_curve(y1, pred_probabilities)\n",
    "    # compute roc auc\n",
    "    #roc_auc = roc_auc_score(y1, pred_probabilities, average = 'micro')\n",
    "    \n",
    "\n",
    "    # compute calibration curve\n",
    "    y, x = calibration_curve(y_test, x1, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array = bin_total(y, x1 , n_bins=10)\n",
    "\n",
    "\n",
    "    print(\"plot curves and save in one png file\")\n",
    "    #save three plots in one png file\n",
    "    fig, (ax1) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    # plot calibration curve LR\n",
    "    ax1.plot(x, y, 'C2',marker='o', linewidth=1, label='LogR')\n",
    "    line = mlines.Line2D([0, 1], [0, 1], color='black')\n",
    "    transform = ax1.transAxes\n",
    "    line.set_transform(transform)\n",
    "    ax1.add_line(line)\n",
    "    ax1.title.set_text(str(plot_name))\n",
    "    ax1.set(xlabel= 'Predicted probability', ylabel= 'True probability in each bin')\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "\n",
    "    # plot calibration curve XGBoost\n",
    "    y, x2 = calibration_curve(y_test, x2, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array2 = bin_total(y, x2 , n_bins=10)\n",
    "    ax1.plot(x2, y, 'C3',marker='o', linewidth=1, label='XGBoost')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    # plot calibration curve RF\n",
    "    y, x3 = calibration_curve(y_test, x3, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array3 = bin_total(y, x3 , n_bins=10)\n",
    "    ax1.plot(x3, y, 'C4',marker='o', linewidth=1, label='RF')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    # plot calibration curve LSTM\n",
    "    y, x4 = calibration_curve(y_test, x4, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array4 = bin_total(y, x4 , n_bins=10)\n",
    "    ax1.plot(x4, y, 'C5',marker='o', linewidth=1, label='LSTM')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    #plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n",
    "y_LR = y_test\n",
    "y_XGB = y_test\n",
    "y_RF = y_test\n",
    "with open('test.npy', 'rb') as f:\n",
    "    #y_lstm = np.load(f)\n",
    "    lstm_prob = np.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibration_together(y_test,LR1_probs,GBT_probs,RF_probs,lstm_prob,  \"Calibration curves for LogR, XGBoost, RF and LSTM\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dist_together (x1,x2, x3,x4, plot_name):\n",
    "\n",
    "    fig, ( ax2) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    # working\n",
    "    ax2.title.set_text(str(plot_name))\n",
    "    ax2.hist(x1, color = \"green\",range=(0, 1), bins=10,label=\"LogR\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x2, color = \"red\", range=(0, 1), bins= 10, label=\"XGBoost\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x3, color = \"purple\", range=(0, 1), bins= 10, label=\"RF\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x4, color = \"brown\", range=(0, 1), bins= 10, label=\"LSTM\", histtype=\"step\", lw=2)\n",
    "    ax2.set(xlabel= 'Mean predicted value', ylabel= 'Count')\n",
    "    ax2.legend(loc=\"upper center\")\n",
    "    \n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_together (LR1_probs,GBT_probs,RF_probs,lstm_prob,  \"Probability distributions for LogR, XGBoost, RF amd LSTM\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# for the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calibration_together (y_test,x3,x4,plot_name):\n",
    "  \n",
    "    #save three plots in one png file\n",
    "    fig, (ax1) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    \n",
    "    line = mlines.Line2D([0, 1], [0, 1], color='black')\n",
    "    transform = ax1.transAxes\n",
    "    line.set_transform(transform)\n",
    "    ax1.add_line(line)\n",
    "    ax1.title.set_text(str(plot_name))\n",
    "    ax1.set(xlabel= 'Predicted probability', ylabel= 'True probability in each bin')\n",
    "    #ax1.legend(loc=\"lower right\")\n",
    "\n",
    "    \n",
    "    # plot calibration curve RF\n",
    "    y, x3 = calibration_curve(y_test, x3, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array3 = bin_total(y, x3 , n_bins=10)\n",
    "    ax1.plot(x3, y, 'C4',marker='o', linewidth=1, label='RF')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    # plot calibration curve LSTM\n",
    "    y, x4 = calibration_curve(y_test, x4, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array4 = bin_total(y, x4 , n_bins=10)\n",
    "    ax1.plot(x4, y, 'C5',marker='o', linewidth=1, label='LSTM')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "\n",
    "    \n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n",
    "\n",
    "with open('test.npy', 'rb') as f:\n",
    "    lstm_prob = np.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "calibration_together(y_test,RF_probs,lstm_prob, \"Calibration curves for RF and LSTM\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dist_together (x3,x4, plot_name):\n",
    "\n",
    "    fig, ( ax2) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    # working\n",
    "    ax2.title.set_text(str(plot_name))\n",
    "    ax2.hist(x3, color = \"purple\", range=(0, 1), bins= 10, label=\"RF\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x4, color = \"brown\", range=(0, 1), bins= 10, label=\"LSTM\", histtype=\"step\", lw=2)\n",
    "    ax2.set(xlabel= 'Mean predicted value', ylabel= 'Count')\n",
    "    ax2.legend(loc=\"upper center\")\n",
    "    \n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_together (RF_probs,lstm_prob, \"Probability distributions for RF amd LSTM\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('padded_lstm.npy', 'rb') as f:\n",
    "    c_lstm_labels = np.load(f)\n",
    "    c_lstm_prob = np.load(f)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calibration_together (y_test,x1, x2,x3,x4, y5, x5, plot_name):\n",
    "    \n",
    "    # compute calibration curve\n",
    "    y, x = calibration_curve(y_test, x1, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array = bin_total(y, x1 , n_bins=10)\n",
    "\n",
    "\n",
    "    print(\"plot curves and save in one png file\")\n",
    "    #save three plots in one png file\n",
    "    fig, (ax1) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    # plot calibration curve LR\n",
    "    ax1.plot(x, y, 'C2',marker='o', linewidth=1, label='LogR')\n",
    "    line = mlines.Line2D([0, 1], [0, 1], color='black')\n",
    "    transform = ax1.transAxes\n",
    "    line.set_transform(transform)\n",
    "    ax1.add_line(line)\n",
    "    ax1.title.set_text(str(plot_name))\n",
    "    ax1.set(xlabel= 'Predicted probability', ylabel= 'True probability in each bin')\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "\n",
    "    # plot calibration curve XGBoost\n",
    "    y, x2 = calibration_curve(y_test, x2, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array2 = bin_total(y, x2 , n_bins=10)\n",
    "    ax1.plot(x2, y, 'C3',marker='o', linewidth=1, label='XGBoost')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    # plot calibration curve RF\n",
    "    y, x3 = calibration_curve(y_test, x3, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array3 = bin_total(y, x3 , n_bins=10)\n",
    "    ax1.plot(x3, y, 'C4',marker='o', linewidth=1, label='RF')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    # plot calibration curve LSTM\n",
    "    y, x4 = calibration_curve(y_test, x4, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array4 = bin_total(y, x4 , n_bins=10)\n",
    "    ax1.plot(x4, y, 'C0',marker='o', linewidth=1, label='LSTM before onset')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    \n",
    "    # plot calibration curve LSTM continuous\n",
    "    y5, x5 = calibration_curve(y5, x5, n_bins=10)\n",
    "    #find out which one are the missing bins\n",
    "    bin_array5 = bin_total(y5, x5 , n_bins=10)\n",
    "    ax1.plot(x5, y, 'C6',marker='o', linewidth=1, label='continuous LSTM')\n",
    "    ax1.add_line(line)\n",
    "    ax1.legend(loc=\"lower right\")\n",
    "    \n",
    "    \n",
    "    \n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n",
    "y_LR = y_test\n",
    "y_XGB = y_test\n",
    "y_RF = y_test\n",
    "with open('test.npy', 'rb') as f:\n",
    "    #y_lstm = np.load(f)\n",
    "    lstm_prob = np.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibration_together(y_test,LR1_probs,GBT_probs,RF_probs,lstm_prob, c_lstm_labels,c_lstm_prob, \"Calibration curves for LogR, XGBoost, RF, LSTM before onset and continuous LSTM\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dist_together (x1,x2, x3,x4, x5, plot_name):\n",
    "\n",
    "    fig, ( ax2) = plt.subplots(1, figsize=(7, 7))\n",
    "    fig.subplots_adjust(wspace=0.3, hspace= 0.3)\n",
    "    \n",
    "    # working\n",
    "    ax2.title.set_text(str(plot_name))\n",
    "    ax2.hist(x1, color = \"green\",range=(0, 1), bins=10,label=\"LogR\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x2, color = \"red\", range=(0, 1), bins= 10, label=\"XGBoost\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x3, color = \"purple\", range=(0, 1), bins= 10, label=\"RF\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x4, color = \"brown\", range=(0, 1), bins= 10, label=\"LSTM before onset\", histtype=\"step\", lw=2)\n",
    "    ax2.hist(x5, color = \"black\", range=(0, 1), bins= 10, label=\"continuous LSTM\", histtype=\"step\", lw=2)\n",
    "    ax2.set(xlabel= 'Mean predicted value', ylabel= 'Count')\n",
    "    ax2.legend(loc=\"upper center\")\n",
    "    \n",
    "    plt.savefig(plot_name+\".png\")\n",
    "    plt.show()\n",
    "    \n",
    "dist_together (LR1_probs,GBT_probs,RF_probs,lstm_prob, c_lstm_prob, \"Probability distributions for LogR, XGBoost, RF, LSTM before onset and continuous LSTM\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Individual graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distribution(LR1_probs, ' LogR no imputation')\n",
    "plt.savefig('dist_LR_NOimp.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classifier_name = \"None vs. Any AKI\"    ###change every time #Moderate vs. Severe #None vs. Any #Others vs. Severe\n",
    "plot_name = \"LogR no imputation\"\n",
    "build_graphs(y_test, LR1_probs, classifier_name, plot_name, \"LR\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XGB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distribution(GBT_probs, ' XGBoost no imputation')\n",
    "plt.savefig('dist_XGB_NOimp.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classifier_name = \"None vs. Any AKI\"    ###change every time #Moderate vs. Severe #None vs. Any #Others vs. Severe\n",
    "plot_name = \"XGBoost no imputation\"\n",
    "build_graphs(y_test, GBT_probs, classifier_name, plot_name, \"XGBoost\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distribution(RF_probs, ' RF no imputation')\n",
    "plt.savefig('dist_RF_NOimp.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "classifier_name = \"None vs. Any AKI\"    ###change every time #Moderate vs. Severe #None vs. Any #Others vs. Severe\n",
    "plot_name = \"RF no imputation\"\n",
    "build_graphs(y_test, RF_probs, classifier_name, plot_name, \"RF\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hyperparameters tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LogR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first default \n",
    "print(datetime.now())\n",
    "# Train binary-classification model with logistic regression\n",
    "print(\"build LR\")\n",
    "model_LR1 = LogisticRegression(penalty = \"l1\", C = 1,fit_intercept = True, solver = 'saga', max_iter= 1000, \n",
    "                                  class_weight= 'balanced', random_state=42)\n",
    "print(\"train LR\")\n",
    "model_LR1.fit(X_train, y_train) \n",
    "# Make predictions\n",
    "print(\"predictions before tuning\")\n",
    "LR1_pred = model_LR1.predict(X_test)\n",
    "# Predict probabilities on testing points\n",
    "LR1_probs = model_LR1.predict_proba(X_test)[::,1]\n",
    "print(datetime.now())\n",
    "\n",
    "performance(y_test,LR1_probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = [0.001, 0.01, 0.1, 1, 10, 100] \n",
    "INTERCEPT = [True, False]\n",
    "penalty = ['l1', 'l2']\n",
    "\n",
    "# Create the random grid\n",
    "random_grid_LR = {'C': C,'fit_intercept': INTERCEPT, 'penalty':penalty}\n",
    "print(random_grid_LR)\n",
    "# Random search training. Use the random grid to search for best hyperparameters\n",
    "# Random search of parameters, using 5 fold cross validation\n",
    "\n",
    "lr_random = RandomizedSearchCV(estimator = model_LR1, param_distributions = random_grid_LR, n_iter = 25, cv = 5, verbose=2, random_state=42, n_jobs = 1)\n",
    "starttime =datetime.now() \n",
    "lr_random.fit(X_train, y_train)\n",
    "print(\"starttime =\", starttime)\n",
    "print(\"endtime =\", datetime.now())\n",
    "\n",
    "# best parameters  (over 30 hours)\n",
    "#'fit_intercept': True, 'C': 0.1 , intercept True, penalty l1.\n",
    "lr_random.best_params_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XGB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counter=collections.Counter(y_train)\n",
    "SCALE_POS_WEIGHT = round(counter[0]/counter[1],1)\n",
    "print(\"build GBT\")\n",
    "#build gradient boost model\n",
    "xgb1_model = xgb.XGBClassifier(\n",
    " learning_rate = 0.1,\n",
    " n_estimators= 100,# the number of boosting rounds or trees to build, equal to XGBoost API num_boost_round\n",
    " max_depth= 5, #should be between 3-10\n",
    " min_child_weight= 1, # a smaller value is chosen because it is a imbalanced class problem and leaf nodes can have smaller size groups\n",
    " gamma= 0,\n",
    " subsample= 0.8,#typical values range between 0.5-0.9\n",
    " colsample_bytree= 0.8,#typical values range between 0.5-0.9\n",
    " objective= 'binary:logistic',\n",
    " nthread= 4,\n",
    " scale_pos_weight= SCALE_POS_WEIGHT,\n",
    " seed= 42)\n",
    "\n",
    "#AL: reg_lambda=1 # default\n",
    "\n",
    "print(\"train GBT\")\n",
    "# train xgb\n",
    "xgb1_model.fit(X_train, y_train, \n",
    "                early_stopping_rounds= 40, \n",
    "                eval_metric= [\"auc\", \"error\", \"logloss\"] , \n",
    "                eval_set= [(X_test, y_test)],\n",
    "                verbose= False)\n",
    "print(\"predictions before tuning\")\n",
    "#Make predictions \n",
    "GBT1_pred = xgb1_model.predict(X_test)\n",
    "GBT1_probs = xgb1_model.predict_proba(X_test)[::,1]\n",
    "performance(y_test,GBT1_probs)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "learning_rate = [0.01,0.1]\n",
    "n_estimators = [20,30,40,50,60,70,80,100]\n",
    "max_depth = [5,7,9]\n",
    "min_child_weight = [1,3,5,7]\n",
    "gamma =  [0, 0.5, 1]\n",
    "\n",
    "\n",
    "# Create the random grid\n",
    "random_grid_xgb = {'learning_rate': learning_rate,\n",
    "               'n_estimators': n_estimators,\n",
    "               'max_depth': max_depth,\n",
    "               'min_child_weight': min_child_weight,\n",
    "               'gamma': gamma}\n",
    "print(random_grid_xgb)\n",
    "\n",
    "# Random search training\n",
    "xgb_random = RandomizedSearchCV(estimator = xgb_model, param_distributions = random_grid_xgb, n_iter = 50, cv = 5, verbose=2, random_state=42, n_jobs = 1)\n",
    "\n",
    "\n",
    "starttime =datetime.now() \n",
    "xgb_random.fit(X_train, y_train)\n",
    "print(\"starttime =\", starttime)\n",
    "print(\"endtime =\", datetime.now())\n",
    "\n",
    "#xgb_random.best_params_\n",
    "\n",
    "# best parameters for full features data set \n",
    "# full features\n",
    "'''\n",
    "{'n_estimators': 70,\n",
    "'min_child_weight': 3,\n",
    "'max_depth': 7,\n",
    "'learning_rate': 0.1,\n",
    "'gamma': 0.5}\n",
    "\n",
    "'''\n",
    "\n",
    "xgb_random.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# search how to do a CV=5 in each model (sklearn function maybe)\n",
    "\n",
    "print(\"build RF\")\n",
    "# building RF (default parameters) to compare to later\n",
    "rf1_model = RandomForestClassifier(n_estimators=100,criterion='gini',max_depth=None,min_samples_split=2,min_samples_leaf=1,\n",
    "                                  min_weight_fraction_leaf=0.0,max_features='auto',max_leaf_nodes=None,min_impurity_decrease=0.0,\n",
    "                                  min_impurity_split=None,bootstrap=True,oob_score=False,n_jobs=None,random_state=None,verbose=0,\n",
    "                                  warm_start=False,class_weight=None,ccp_alpha=0.0,max_samples=None)\n",
    "\n",
    "\n",
    "print(\"train RF\")\n",
    "# train RF\n",
    "rf1_model.fit(X_train, y_train)\n",
    "# to see how the forest looks like with default parameters\n",
    "\n",
    "# depth of trees is around 44-93 for the full set\n",
    "#[estimator.tree_.max_depth for estimator in rf_model.estimators_]\n",
    "\n",
    "print(\"predictions\")\n",
    "#Make predictions \n",
    "RF1_pred = rf1_model.predict(X_test)\n",
    "RF1_probs = rf1_model.predict_proba(X_test)[::,1] # probability of class 1\n",
    "\n",
    "performance(y_test,RF1_probs)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Number of trees in random forest\n",
    "n_estimators = [int(x) for x in np.linspace(start = 100, stop = 1000, num = 10)]\n",
    "# Number of features to consider at every split\n",
    "max_features = ['auto']\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [int(x) for x in np.linspace(10, 60, num = 10)]\n",
    "max_depth.append(None)\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [2, 5, 10]\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [2, 4]\n",
    "# Method of selecting samples for training each tree\n",
    "bootstrap = [True, False]\n",
    "# Create the random grid\n",
    "random_grid = {'n_estimators': n_estimators,\n",
    "               'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf,\n",
    "               'bootstrap': bootstrap}\n",
    "print(random_grid)\n",
    "\n",
    "# Random search training\n",
    "# Use the random grid to search for best hyperparameters\n",
    "# Random search of parameters, using 3 fold cross validation\n",
    "rf_random = RandomizedSearchCV(estimator = rf_model, param_distributions = random_grid, n_iter = 100, cv = 5, verbose=2, random_state=42, n_jobs = 1)\n",
    "\n",
    "\n",
    "#Fit the random search model  \n",
    "starttime =datetime.now() \n",
    "rf_random.fit(X_train, y_train)\n",
    "print(\"starttime =\", starttime)\n",
    "print(\"endtime =\", datetime.now())\n",
    "\n",
    "rf_random.best_params_\n",
    "\n",
    "# max optimized\n",
    "'''\n",
    "{'n_estimators': 700,\n",
    "'min_samples_split': 10,\n",
    "'min_samples_leaf': 2,\n",
    "'max_features': 'auto',\n",
    "'max_depth': 40,\n",
    "'bootstrap': False}\n",
    "'''\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
